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Abstract 

The physics of crystalhne membranes, i.e. fixed-connectivity surfaces em- 
bedded in three dimensions and with an extrinsic curvature term, is very rich 
and of great theoretical interest. To understand their behavior, numerical simu- 
lations are commonly used. Unfortunately, traditional Monte Carlo algorithms 
suffer from very long auto-correlations and critical slowing down in the more 
interesting phases of the model. In this paper we study the performance of 
improved Monte Carlo algorithms for simulating crystalline membrane, such as 
hybrid overrelaxation and unigrid methods, and compare their performance to 
the more traditional Metropolis algorithm. We find that although the overre- 
laxation algorithm does not reduce the critical slowing down, it gives an overall 
gain of a factor 15 over the Metropolis algorithm. The unigrid algorithm does, 
on the other hand, reduce the critical slowing down exponent to z sa 1.7. 



1 Introduction 

When using Monte Carlo methods to study physical systems one is usually faced 
with the problem of critical slowing- down (CSD) in the critical region were a typ- 
ical correlation length of the system diverges. That is, the auto-correlation time 
of traditional Monte Carlo algorithms, the time it takes to generate "statistically 
independent" configurations, grows rapidly with system size. For example, for the 
well known Metropolis and heat-bath algorithms this time grows linearly with the 
system size, making simulations of large systems prohibitively time consuming. 

It is thus very important to construct new Monte Carlo algorithms that can 
reduce CSD. Much progress has been made in this direction; examples of improved 
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algorithms are, to name a few: Adler's overrelaxation Fourier acceleration |2|, 
multigrid and cluster algorithms 0] . Those algorithms have been applied success- 
fully to variety of models and, in some special cases, have eliminated CSD altogether. 
This, though, is usually accomplished only for relatively simple models, such as free 
field theories or spin models with simple interactions. For more complicated models, 
where sophisticated methods are harder to implement, the improvement is usually 
somewhat less. 

In this paper we study the performance of two such improved algorithms, 
overrelaxation and multigrid, in simulations of crystalline membranes. Crystalline 
membranes are internally rigid surfaces, embedded in three dimensions, with an ex- 
trinsic curvature term. They exhibit a phase transition between a high-temperature 
crumpled and a low-temperature flat phase. It is especially in the flat phase, and at 
the crumpling transition, that Monte Carlo simulations with traditional algorithms 
suffer from very long auto-correlations. 

As with many interesting models, the Hamiltonian of a crystalline membrane 
is too complicated for a direct implementation of the methods we want to employ. 
Some simplifications have to be made. For the overrelaxation we use a quadratic 
approximation to the Hamiltonian when choosing a new trial position — this requires 
an additional accept-r eject step to restore detail balance. This is usually referred to 
as hybrid overrelaxation. Instead of multigrid we use a simpler implementation, the 
unigrid algorithm, in which a coarsening transformation of the field configuration is 
not needed — the fields at the original (fine) lattice are simulated at all levels. 

In addition to the performance of the algorithms, we also examine the impor- 
tance of randomness in the updating procedure, i.e. the order in which the fields 
are updated. According to standard folklore, too much randomness in the updating 
procedure increases CSD as the system takes a drunkard's walk through the phase 
space. Too little noise, on the other hand, increases CSD as well since the system is 
too weakly ergodic. It is thus important to tune the amount of randomness in the 
algorithm appropriately. 

The paper is organized as follows. In Section 2 we describe the particular 
model of a crystalline membrane we study, discuss the problems of the simulations 
and the performance of the Metropolis algorithm. In Section 3 we describe the 
hybrid overrelaxation algorithm and our approximation. In Section 4 we test the 
performance of the unigrid algorithm. Finally, in the Section 5, we compare the 
overall performance and merits of those different methods and comment on possible 
further improvements and applications. 

2 A model of crystalline membranes 

The model we simulate is a simple discretization of a phantom (non self-avoiding) 
crystalline membrane, inspired by the Polyakov action for Euclidean strings with 
extrinsic curvature A discrete crystalline membrane is described by a reg- 
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ular two-dimensional triangulation embedded in three-dimensional space where it 
is allowed to fluctuate. The Hamiltonian is composed of two terms: a pair poten- 
tial between neighboring nodes and a bending energy term. As a pair potential 
we use a simple Gaussian spring potential, and we model the bending energy by a 
ferromagnetic interaction between neighboring normals to the faces of the surface: 

= + n^il - fia ■ fib). (1) 

(ij) {ab) 

Here i,j label the intrinsic position of the nodes, fi is the corresponding position in 
the embedding space, is a normal to a triangle a, and k is the bending rigidity. 
The partition function is given by the trace of the Boltzmann weight over all possible 
configurations of the embedding variables {r}: 

Z = J [df ] 5(fe„) exp {-n[r] ) . (2) 

The center of mass of the membrane, fcm, is kept fixed to eliminate the transla- 
tional zero mode. As there is no self-avoidance term in the Hamiltonian, this model 
describes a phantom surface. 

This model has been studied extensively with numerical methods 1^, 0, |8|, |^. 
It has been found to have a high-temperature crumpled (disordered) phase and a 
low-temperature flat (ordered) phase, separated by a continuous phase transition — 
the crumpling transition. The behavior of the system in the flat phase is governed 
by an infrared stable fixed point at «; = oo; the whole flat phase is critical. The 
existence of an ordered phase in a two-dimensional system with a continuous symme- 
try and short range interactions is remarkable, given the Mermin- Wagner theorem. 
What stabilizes the flat phase are the out-of-plane fluctuation of the membrane 
that couple to the in-plane phonon degrees of freedom due to non-vanishing elas- 
tic moduli. Bending of the membrane is necessarily accompanied by an internal 
stretching. By integrating out the phonon degrees of freedom, one is left with an 
effective Hamiltonian with long-range interactions between the Gaussian curvature 
fluctuations. 

Most numerical simulations so far have used either local updating methods, 
usually the Metropolis algorithm, or Fourier acceleration. The Metropolis algorithm, 
apart from suffering from CSD, has very long auto-correlations both in the flat phase 
and close to the crumpling transition. To establish this we have simulated the model 
Eq. (|2|), using the Metropolis algorithm, on a L x L square lattice, L ranging from 
8 to 64, and with periodic boundary conditions. We choose to simulate the model 
in the flat phase, at k = 1.1, were we know from previous simulations that the 
auto-correlation time is indeed very long Q. 

To estimate the auto-correlations we measure the square radius of gyration: 
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Table 1: The integrated auto-correlation time r (in number of sweeps), together with the 
CPU-time per sweep (in ms). Column (a) is for a lexicographic update of nodes, while 
(b) is for random updating. From a linear fit to Eqs. (|^) and (|^) we get the exponents Za 
and Zg and the corresponding amplitudes Aa and As- 





(a) 




(b) 




L 


T 


Ts 


T 


Ts 


8 


219(15) 


1.112 


227(12) 


1.036 


12 


546(35) 


2.784 


567(30) 


2.610 


16 


1153(90) 


4.900 


1123(80) 


4.688 


24 


2714(150) 


10.86 


2534(220) 


11.01 


32 


4049(210) 


22.23 


4527(260) 


21.68 


48 


11443(600) 


52.33 


10347(550) 


50.60 


64 


20500(1200) 


97.59 


19900(1300) 


98.18 


z 


2T61(34) 


2.143(23) 


2.128(34) 


2.188(31) 


A 


2.57(29) 


0.0129(10) 


2.83(28) 


0.0111(13) 



i.e. the linear extent of the membrane in the embedding space. This is usually 
the "slowest mode" of the system. From this we construct the normalized auto- 
correlation function 

,(.) = w+^mt))-^(^.)\ 

and the integrated auto-correlation time (measured in units of sweeps) 

^ s=l 

This is shown in Table 1. The errors on the auto-correlation times are estimated 
from 10 to 20 independent runs, each few hundred r long. 

As mentioned in the introduction, we also want to understand the effect of 
randomness in the updating on performance. Thus we have repeated the simulations 
for two different updating schemes: the nodes are sampled either at random or they 
are traversed in a lexicographic order. Lexicographic order means that the node 
at (intrinsic) position x is always updated before x + Ci {i = 1,2), except at the 
boundaries. For a free field theory, updated with the overrelaxation algorithm, only 
the latter scheme reduces CSD [p^]. 

As an estimate of CSD we define the dynamical critical exponent Za using 
finite size scaling: 

r « (6) 
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Here is some characteristic length scale of the system; in our simulations, where 
the model is critical, this is the intrinsic linear extent L. 

To evaluate the performance of the algorithm, we also have to take into account 
the amplitude Aa and how the computational work performed per sweep, measured 
in CPU-time Tg, scales with system size. This we include in Table 1. Similarly to 
Eq. (^) we define an exponent Zg-. 

Ts ^ Ase% (7) 

where t is now measured in "real" time (in ms). Combining these exponents and 
amplitudes, the performance — the total "cost" of the algorithm — is given by 
T = AgAat^"'^^" ■ For the Metropolis algorithm we get: 

_ f 0.0332(45) L^-304(57) lexicographic updates, , . 

~ \ 0.0314(48) l4.316(67) random updates. ^ ' 

For this algorithm, the order in which nodes are updated is irrelevant. 



3 Hybrid over relaxation 

Overrelaxation was introduced as a generalization of the heat-bath algorithm for 
models with multi-quadratic actions Q. In the original formulation the new value 
of a field 4)i is chosen to be negatively correlated with the old value. Given a multi- 
quadratic action, 

S = LJ {(pi — J^i[(f)j^i])^ + { terms independent oi (pi }, (9) 
one chooses a local update of the field (pi as 



cp, ^ (P', = {I- c)(Pi + + ^l^KJl^ . (10) 

1^ is a Gaussian random variable of unit variance and C is a relaxation parameter. 
This update fulfills detail balance for < C < 2; for ^ = 1 it reduces to the standard 
heat-bath algorithm, while for C = 2 the field evolution becomes deterministic and 
conserves energy (a micro-canonical simulation). In the latter case, in order to 
restore ergodicity, some amount of standard ergodic updates have to be included. 

This method has been applied successfully to variety of models; its success 
based on it suppressing the usual random walk behavior of local updating algo- 
rithms |]l^]. In order to achieve the greatest reduction of CSD, both the relaxation 



parameter and the noise in the updating procedure, should be fine-tuned ||10|, |14|. 
Unfortunately, the usefulness of the method has been limited by its restriction to 
multi-quadratic systems. 
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A number of generalizations of the overrelaxation have been proposed |12|. 
They usuahy involve a transformation of the Boltzmann distribution to the appro- 
priate form Eq. (9), and the introduction of an accept-reject step to ensure detail 
balance. This has the disadvantage that the accept-reject step can enhance random 
walk behavior by the algorithm and, in addition, the rejection probability usually 
depends on some characteristics of the model and may not be adjustable to a rea- 
sonable value. This is nevertheless the approach we will use. 

We make a quadratic approximation to the Hamiltonian Eq. ([l|) and then 
apply hybrid overrelaxation (with (" = 2). We treat the non-linear bending energy 
term by assuming that normalization in the denominator is constant for all the 
triangles, i.e. we write the normals as 

-. ^-.^ ^ {rj -fj)x {rj - ffc) ^ fj xfk-fjX (fj + ffc) , ^^^^ 

were {i,j, k} are the nodes defining triangle a. Then, a quadratic approximation to 
the Hamiltonian can be written as 

= X! ~ '^jf ~ [^i X n-fix {fj + fk)] ■ [n xfi-fix (ffc + fi)] , 

(ij) (ab) 

(12) 

where the triangles a = k) and b = {i, k, I) are adjacent. Since the approximate 
Hamiltonian is quadratic in rj, we can write: 

Ti.A = fi • (Mfi) + C ■ fi + {terms independent of fj}. (13) 

The matrix M and the vector C are easily computed: 



«-^i:E('-f-fO(o?.--f 

j=lb^a 



and 



C^"^ = - E {rf {r, • (2r-;-i + 2f,+i - - q,.,) - (f,+i + f,_i) 



+ i2^^^(-l^^-rl?\ . (15) 



The indices a and b label the component of the fields in the embedding space, and 
the index j labels the neighbors of node i, including its next-to-nearest neighbors 
Ijj, in a cyclic manner. 
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Figure 1: The integrated auto-correlation time t, vs. the "normalization" parameter A, 
for the hybrid overrelaxation algorithm. This is shown both for a random (squares) and a 
lexicographic (circles) updating. The lattice size is 16^. The dashed line is the corresponding 
acceptance rate p in the Metropolis test. 



The constant energy surface, TCAi^i) = k, is a multi-quadratic function — 
in our case an ellipsoid in the embedding space. To find the new (overrelaxed) 
embedding position r^', we can diagonalize the matrix M and apply overrelaxation 
to the transformed variables. This involves some amount of calculation; a quicker 
and STifRcient method is to apply overrelaxation in a random sequence to each of 
the embedding positions r|"\ a = 1,2,3. Once the trial position has been chosen, 
it is accepted or rejected according to a Metropolis test. 

An important feature of this approximation is the parameter A. Although 
introduced as a substitution for the normalization of the normals, it can be tuned to 
optimize the performance of the algorithm by minimizing the rejection probability 
in the Metropolis test. 

We have applied this method in the fiat phase {k = 1.1) for both random and 
lexicographic updating. To ensure ergodicity we also include a random amount of 
standard Metropolis updates (about 20%). In Fig. 1 we show the integrated auto- 
correlation time vs. A, for a lattice size L = 16. We also include the corresponding 
acceptance rate in the Metropolis test. In this case, contrary to the Metropolis 
algorithm, lexicographic updating reduces r by about 30%, independent of system 
size. More important, for a suitable choice of A, r is reduced by a factor of 15 
relative to the Metropolis algorithm. As it might be expected, the optimal choice of 
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Table 2: Same as Table 1, except the algorithm used is the approximate hybrid overrelax- 
ation. Again (a) corresponds to lexicographic and (b) to random updating. The optimal 
value of the parameter A is also included. 







(a) 






(h) 




L 


A 


T 


T 

^ sweep 


A 


T 


sweep 


8 


1.08 


18.5(7) 


1.247 


1.08 


31.9(1.2) 


1.136 


12 


1.23 


36.6(1.1) 


3.078 


1.22 


66.9(2.0) 


2.766 


16 


1.38 


70.6(8.0) 


5.735 


1.35 


113.3(6.0) 


5.122 


24 


1.57 


150(11) 


13.42 


1.58 


279(18) 


11.83 


32 


1.92 


269(25) 


27.27 


1.94 


484(42) 


24.12 


48 


2.68 


640(40) 


62.69 


2.60 


1096(90) 


53.19 


64 


3.09 


1120(105) 


118.7 


3.20 


2600(180) 


103.5 


z 




1.990(31) 


2.190(14) 




2.065(30) 


2.163(18) 


A 




0.275(24) 


0.0132(6) 




0.405(34) 


0.0127(8) 



A corresponds to maximizing acceptance rate in the Metropolis test. Surprisingly 
this value is much higher than one would expect from the average length of the 
un-normalized normals {{\nun\) ~ 0.3 for k = 1.1). This implies that, in order to 
enhance the acceptance rate, it is convenient to suppress the bending energy term 
in the approximation. 

We have repeated this analysis for lattices sizes L = 8 to 64. In Table 2 we 
show the optimal values of A, the corresponding integrated auto-correlation time 
and the CPU-time for a sweep. Prom this data we extract, as before, the exponents 
Za and Zs and the amplitudes Aa and As-, and obtain the following performance: 

_ f 0.00364(36) L'^-i^^^^^) lexicographic updates, , . 

^ ~ I 0.00514(30) l4-228{48) random updates. ^ ' 

Although hybrid overrelaxation does not reduce CSD, it gives an improvement of one 
order of magnitude over the Metropolis algorithm, provided the nodes are updated 
in a lexicographic order and the "normalization" parameter A is properly adjusted. 

4 A unigrid Monte Carlo algorithm 

The critical slowing down of traditional Monte Carlo algorithms arises mainly from 
the fact that the update is local, and thus the system takes a random walk through 
the configuration space. This can be improved by using collective mode (non-local) 
updating such as multigrid methods fs], 13|. The basic idea is to consider a sequence 



of coarser lattices (levels) in addition to the original lattice. At each level the system 
is updated using traditional methods but, as this is repeated recursively at all length 
scales, the long wavelength modes are equilibrated much faster. 
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There are several basic ingredients to a multigrid algorithm: a restriction 
operator and the corresponding interpolation operator, or kernel, are needed to 
map the system onto the coarser lattices and back; an updating algorithm, such as 
Metropolis, is applied at each level; finally, one has to choose how to traverse the 
different levels. 

For a crystalline membrane it is not possible to construct an exact interpo- 
lation operator between different levels due to the complexity of the Hamiltonian 
Eq. (H). This problem can be circumvented by using an alternative implementation, 
the unigrid method, in which the coarse lattices are simply defined as subdivisions 
of the original one; the update of a coarse lattice acts on blocks of the original fields. 
For the update, the choice is usually between a piecewise constant or a piecewise 
linear kernel. A piecewise kernel simply shifts all of the fields in the block by a uni- 
form value. A piecewise linear kernel shifts the fields by a value linearly interpolated 
between zero, at the boundary, and a maximum value at the center of the block. 
The shift operation is one of the global symmetries of the system. 

Several considerations should be made in choosing a kernel. The piecewise 
linear kernel has the advantage that the acceptance rate of the proposed moves does 
not depend on the block size. For a crystalline membrane, on the other hand, this is 
outweighed by the computational cost which, as all the normal-normal interactions 
in the block have to be recalculated, grows linearly with the block size. Hence a 
piecewise constant kernel is preferable, since the sole contribution to the energy 
change comes from the boundary. 

We parameterize a non-local change of the configuration, when we shift a 
block A/j at level k, as: 

^/ ^ i ri + ekx if ie Ak .^^^ 
* I fi otherwise, 

where x is some normalized random noise and is the amplitude of the shift. A 
necessary prerequisite for the unigrid method to be efficient is that the energy cost 
does not grow too fast with the perimeter of the block. Stated differently: in 
order to maintain a constant acceptance rate in the Metropolis test, the amplitudes 
have to be scaled like 

~ ^r- (18) 

If the exponent a is too big, i.e. of order unity, it is unlikely that any multigrid 
algorithm will reduce CSD. 



Following the analysis of [15 1, we can estimate whether the unigrid update has 
a chance of improving the dynamical behavior in the case of a crystalline membranes. 
Assuming that the probability distribution is nearly Gaussian, one can approximate 
the acceptance rate by 

J^(e) ~ erfc(\/^/2), (19) 

where h = {ATC) is the average change in energy and e is the amplitude of the 
proposed move. For a piecewise constant kernel we take the Hamiltonian Eq. (1), 
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Table 3: The auto-correlation and CPU-times for the unigrid algorithm. Results are shown 
both for V and IT-cycles. 





V -cycle 


W -cycle 


L 


T 


Ts 


T 


Ts 


8 


17.1(0.4) 


3.123 


14.0(0.3) 


4.195 


16 


45.8(3.2) 


16.325 


30.1(0.9) 


25.904 


32 


115.2(6.1) 


107.18 


51.0(2.1) 


203.24 


64 


269.3(19.1) 


524.70 


96.8(8.7) 


1324,25 


z 


1.349(28) 


2.489(55) 


0.955(27) 


2.788(45) 


A 


1.040(77) 


0.0175(31) 


3.96(14) 


0.0123(18) 



insert Eq. ([T^) and expand the change in the energy in powers of e: 

(AW({r-;};x)) = e(Fi({r-;};x)) + {F2{{n};x)) + 0{e^). (20) 

The key observation is that under a global sign change, fi — > —fi Vz, the func- 
tion Fi changes sign (as each term is a product of odd number of fields r^), hence 
(-Fi({rj}; x)) = 0. The leading contribution to h is therefore proportional to e^. At 
the same time, the number of terms contributing to Eq. ( pO| ) depends linearly on 
the length of the boundary. Therefore, in order to maintain a constant acceptance 

— 1/2 

rate, one should scale the amplitudes as e ~ . This agrees with our numerical 

simulations where we find a = 0.52(1). 

Another free parameter in the unigrid method is the relative frequency with 
which different levels are updated. One must strike a balance between the effec- 
tiveness of block moves and their relative computational cost. Two general schemes 
are used in the literature; the V and T^-cycles. In a F-cycle the levels are up- 
dated consecutively, from the finest to the coarsest and back, whereas a M^-cycle 
recursively applies a V-cyc\e at each visited level, spending more time in updating 
coarser levels. For multigrid algorithms, where the lattices size decreases between 
levels, a VK-cycle is preferable, provided that the interpolations between levels is 
not too time consuming. For a unigrid algorithm the computational cost increases 
with the block size, as discussed above, and, depending on the exponent a, a F or 
W— cycle will be advantageous. For a piecewise constant kernel the computational 
cost scales like Lb and a W -cycle might be advantageous, while for a piecewise 
linear kernel the computational cost scales linearly with the area of the block and a 
F-cycle would be better. 

As before, we have simulated a crystalline membrane in the flat phase and 
for lattice sizes L = 8, 16, 32 and 64, using the unigrid algorithm and updating the 
system at each level with the Metropolis algorithm. We repeated the simulation 
both for a V and T^-cycle. In Table 3 the show the measured value of the auto- 
correlation and CPU-times, from which we determine the overall performance of the 
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algorithm: 

_ / 0.0182(35) l3-838(84) y_cycle, , . 

- \ 0.0242(39) l3-743(72) vF-cycle. ^ ^ 

In both cases the unigrid algorithm reduces CSD, albeit not very much, but enough 
to make it worthwhile for large membranes. For the VF-cycle this implies a dynam- 
ical CSD exponent z ^ 1.7, although this value is probably to large, as we observe 
strong finite size effects in the fit to Eq. (^); if we exclude the smallest volume 
(L = 8), we get z ~ 1.6. In conclusion, although the VF-cycle is more time consum- 
ing per sweep, its performance is better than l^-cycle in simulations of membranes 
of size L > 20. 



5 Discussion 



Comparing the performance of these different algorithms, Eqs. (^), (|T6|), and (^T]), 
we see that, in the simulation of crystalline membranes on realistic lattice sizes 
(L = 10 to 200), both the hybrid overrelaxation and the unigrid algorithm reduce 
the computational cost by an order of magnitude over the Metropolis algorithm. 
As the unigrid algorithm also reduces the dynamical exponent z, especially using a 
-cycle, it is clearly the best choice for large membranes. In the particular case we 
studied in this paper, large means L > 50, although this value may depend on the 
simulation parameters (i.e. k). 

We would also like to emphasize that, in order to achive optimal performance 
of the hybrid overrelaxation, it is imperitive to adjust the noise in the updating pro- 
cedure (to use lexicographic updating), and to tune the "normalization" parameter 
A appropriatly. 

An alternative algorithm used to simulate crystalline membranes is a com- 
bination of Langevin updates with Fourier acceleration ||^, |^. This algorithm is 
known to substantially reduce CSD, although the gain is lost to some extent in the 
large computational overhead. This method is also complicated by systematic errors 
induced by using a finite time step At; this necessitates an extrapoltion to At = 0, 
which can itself become time consuming. Nevertheless, it would be interesting to 
know how well this algorithm performes, in realistic simulations, compared to the 
algorithms we have studied in this paper. Unfortunately, we do not have an esti- 
mate of its performance in similar conditions (e.g., using the same computers) for 
comparison. 

An obvious extension of the methods studied in this paper, is to combine 
hybrid overrelaxation with the unigrid algorithm. It is possible to maximize the shift 
of a block in a unigrid update, by choosing it in a deterministic and energy preserving 
way, improving the performance even further. This is though more complicated to 
implement and it is under investigation. 

There are few application in which the hybrid overrelaxation might be more 
advantageous compared to the unigrid algorithm. Overrelaxation can be parallelized 
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in a straightforward manner, although it could be difficult to define a lexicographic 
update in that case. It also easy to adapt hybrid overrelaxation to modified versions 
of the model like, for example, self-avoiding crystalline membranes, which are of 
great physical interest. In that proposed updated is rejected if it leads to 

self-intersection of the membrane. Intuitively, a non-local change, like one proposed 
by the unigrid algorithm, is more likely to be rejected — hybrid overrelaxation might 
turn out to be more effective. 

Finally, we would like to point out that both these updating algorithms can be 
used in simulations of fluid membranes with extrinsic curvature |16|. In that case, 
the surface fluctuates in the embedding space, and its connectivity matrix changes 
dynamically. It is known that simulations of fluid membranes suffer from even 
longer auto-correlation times than their crystalline counterparts. For fluid mem- 
branes, the random nature of the lattice causes some complications in implementing 
hybrid overrelaxation and unigrid algorithms. For example, it is not obvious how 
to define lexicographic ordering. A possible method would be to propagate the 
updates outwards from a randomly chosen node, i.e. by traversing the lattice in 
steps of increasing geodesic distance from a marked point. For the unigrid algo- 
rithm the problem is that a random lattice cannot be divided into regular blocks. 
Again, blocks could be defined as nodes within a given geodesic distance from some 
randomly chosen node. 
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